﻿#region Copyright information
// 
// Copyright © 2005-2013 Yongkee Cho. All rights reserved.
// 
// This code is a part of the Biological Object Library and governed under the terms of the
// GNU Lesser General  Public License (LGPL) version 2.1 which accompanies this distribution.
// For more information on the LGPL, please visit http://bol.codeplex.com/license.
// 
// - Filename: PoissonDistribution.cs
// - Author: Yongkee Cho
// - Email: yongkeecho@gmail.com
// - Date Created: 2012-09-06 11:39 AM
// - Last Modified: 2013-01-25 3:59 PM
// 
#endregion
using System;
using System.Collections.Generic;
using System.Linq;
using BOL.Maths.Functions;
using BOL.Algorithms.Optimization;

namespace BOL.Maths.Distributions
{
    /// <summary>
    /// Represents a Poisson distribution with mean mu.
    /// <see>Mathematica: ref/PoissonDistribution</see>
    /// <see>http://en.wikipedia.org/wiki/Poisson_distribution</see>
    /// </summary>
    public sealed class PoissonDistribution : IUnivariateDistribution<int>, IEquatable<PoissonDistribution>
    {
        #region Private variables

        private double _lambda;

        #endregion

        #region Public properties

        /// <summary>Gets or sets the parameter, lambda of a Poisson distribution.</summary>
        public double Lambda { get { return _lambda; } set { if (value <= 0) throw new ArgumentOutOfRangeException("value", "lambda is expected to be a positive real number."); _lambda = value; } }

        /// <summary>Gets the range of a Poisson distribution.</summary>
        public IRange<int> Domain { get { return new Range<int>(0, int.MaxValue); } }

        /// <summary>Gets the mean of a Poisson distribution.</summary>
        public double Mean { get { return _lambda; } }

        /// <summary>Gets the median of a Poisson distribution.</summary>
        public double Median { get { return Math.Floor(_lambda + 1 / 3.0 - 0.02 / _lambda); } }

        /// <summary>Gets the mode of a Poisson distribution.</summary>
        public double Mode { get { return Math.Ceiling(_lambda) - 1; } }

        /// <summary>Gets the variance of a Poisson distribution.</summary>
        public double Variance { get { return _lambda; } }

        /// <summary>Gets the skewness of a Poisson distribution.</summary>
        public double Skewness { get { return 1 / Math.Sqrt(_lambda); } }

        /// <summary>Gets the kurtosis of a Poisson distribution.</summary>
        public double Kurtosis { get { return 1 / _lambda; } }

        /// <summary>Gets the entropy of a Poisson distribution.</summary>
        public double Entropy { get { throw new NotImplementedException(); } }

        #endregion

        #region Constructor

        /// <summary>
        /// Instantiates a Poisson distribution.
        /// </summary>
        /// <param name="lambda"></param>
        public PoissonDistribution(double lambda)
        {
            if (lambda <= 0)
                throw new ArgumentOutOfRangeException("lambda", "lambda is expected to be positive.");

            _lambda = lambda;
        }

        #endregion

        #region ICloneable implementation

        public PoissonDistribution Clone()
        {
            return new PoissonDistribution(_lambda);
        }

        object ICloneable.Clone()
        {
            return Clone();
        }
        
        #endregion

        #region Public methods

        /// <summary>
        /// Returns the Poisson density given value.
        /// </summary>
        /// <param name="value"></param>
        /// <returns>Poisson density</returns>
        public double Pdf(int value)
        {
            if (value < 0)
                throw new ArgumentOutOfRangeException("value", "value is expected to be a positive integer.");

            return Math.Exp(-_lambda) * Math.Pow (_lambda, value) / CountingFunctions.Factorial(value);
        }

        /// <summary>
        /// Returns the cummulative Poisson density given value.
        /// </summary>
        /// <param name="value"></param>
        /// <returns>cummulative Poisson density</returns>
        public double Cdf(int value)
        {
            return DistributionFunctions.IncompleteGamma((value + 1), _lambda) / CountingFunctions.Factorial(value);
        }

        /// <summary>
        /// Returns the value of the Poisson distribution for given probability.
        /// </summary>
        /// <param name="p">probability</param>
        /// <returns></returns>
        public int Quantile(double p)
        {
            //return Math.Round(IRootFinder.Find(p, _lambda, 0.0, double.MAX_VALUE));

            throw new NotImplementedException();
        }

        /// <summary>
        /// Returns a random value from the Poisson distribution.
        /// <see>http://people.csail.mit.edu/milch/blog/index.html</see>
        /// <see>Anuj Kumar's course page for IEOR E4404 at Columbia University, http://www.columbia.edu/~ak2108/ta/summer2003/poisson1.c</see>
        /// </summary>
        /// <returns>random value</returns>
        public int Sample(Random r)
        {
            var prob = Math.Exp(-_lambda);
            var cumProb = prob;
            var i = 0;
            var u = r.NextDouble();

            while (cumProb < u)
            {
                i++;
                prob *= _lambda / i;
                cumProb += prob;
            }

            return i;
        }

        /// <summary>
        /// Estimates parameters of the Poisson distribution using maximum likelihood.
        /// <see>http://en.wikipedia.org/wiki/Poisson_distribution#Maximum_likelihood</see>
        /// </summary>
        public void MaximumLikelihoodEstimate(IEnumerable<int> source)
        {
            if (source == null)
                throw new ArgumentNullException();

            Lambda = source.Average();
        }

        /// <summary>
        /// Estimates parameters of the Poisson distribution using Bayesian methods.
        /// <see>http://en.wikipedia.org/wiki/Poisson_distribution#Bayesian_inference</see>
        /// </summary>
        public void BayesianEstimate(IEnumerable<int> source, params dynamic[] priorParameters)
        {
            //if (source == null)
            //    throw new ArgumentNullException();

            //var prior = priorParameters.Cast<double>().ToList();
            //var pp = new PoissonPosterior(data, prior[0], prior[1]);
            //var optimizer = new IOptimizerByValue(pp, 5000, 1e-5, true, 0.1);
            //optimizer.Optimize();
            //return new PoissonDistribution(pp.Parameters[0]);

            throw new NotImplementedException();
        }

        #region Poisson posterior function 

        private class PoissonPosterior : ObjectiveFunctionBase
        {
            #region Private variables

            private readonly GammaDistribution _gamma; // a Gamma distribution

            #endregion 

            #region Constructor

            public PoissonPosterior(IEnumerable<double> source, double alpha, double beta)
                : base(1)
            {
                var list = source.ToList();
                _gamma = new GammaDistribution(alpha + list.Average(), beta + list.Count());
            }

            #endregion

            #region Public methods

            public override double Value(double[] parameters)
            {
                return _gamma.Pdf(parameters[0]);
            }
            
            #endregion
        }

        #endregion

        #endregion

        #region IEquatable<PoissonDistribution> implementation

        public bool Equals(PoissonDistribution other)
        {
            return _lambda.Equals(other._lambda);
        }

        #endregion

        #region Object overriden

        public override int GetHashCode()
        {
            return Lambda.GetHashCode();
        }

        public override bool Equals(object other)
        {
            if (other == null)
                throw new ArgumentNullException("other");

            if (!(other is PoissonDistribution))
                throw new InvalidCastException("The 'other' argument is not a PoissonDistribution object.");
            
            return Equals(other as PoissonDistribution);
        }

        public override string ToString()
        {
            return String.Format("Poisson(lambda = {0})", _lambda);
        }

        #endregion
    }
}
